*****************************************************************
* Replication directory for                                   ***
* Prime locations                                             ***
* by Gabriel M. Ahlfeldt, Thilo N.H. Albers, Kristian Behrens ***
* Published in American Economic Review: Insights             ***
*****************************************************************
* 01/2025
* Stata
version 17.0

* This do file generates gradients showing how various outcomes change in 
* distance from prime locations within Global Cities


* Part I: Gradients by urban spatial structure *********************************

* Load basline grid data
	u "$temp/125CitiesPLs/grid125_PL_output.dta", clear
	* drop cities that are not in final sample due to lack of data 	
	foreach num of numlist 18 39 54 113 {
	drop if metro_id == `num'
	}	
									
* classify mono, duo-, and polycentric
	gen N_PLs_city=.
	levelsof metro_id, local(metro_ids)
	foreach x of local metro_ids{
	distinct PLID_metro if metro_id==`x'
	replace N_PLs_city=r(ndistinct) if metro_id==`x'
	}
									 
* Merge distances
	merge 1:1 grid_x grid_y using "$temp/grid125_PL_distances"
	drop if _m == 2
	drop _m
								 
* Merge SNL and Starbucks data
	merge 1:1 grid_x grid_y using "$data_125cities/STARBUCKS_SNL/grid_StarSNL.dta"
	drop if _m == 2
	drop  _m
	distinct metro_id if snl_count!=0 & snl_count!=.
	local N_snl_c=r(ndistinct)
	distinct metro_id if star_count!=0 & star_count!=.
	local N_star_c=r(ndistinct)
									
* MERGE TWITTER DATA
	merge 1:1 grid_x grid_y using "$data_125cities/TWITTER/DATA.dta",  
	tab _m
	drop if _m == 2
	drop _m
	distinct metro_id if twittercount!=0 & twittercount!=.
	local N_tc_c=r(ndistinct)
									
* MERGE EMPORIS DATA 
	merge 1:1 grid_x grid_y using "$data_125cities/EMPORIS/EMPORIS_VALIDATION_maxHeight.dta",  
	tab _m
	drop if _m == 2
	drop _m
	distinct metro_id if twittercount!=0 & twittercount!=.
	local N_tc_c=r(ndistinct)
											
* Merge photo data
	merge 1:1 grid_x grid_y using "$data_125cities/PHOTOS/PHOTOCOUNT_250M.dta" ,   	
	tab _m
	drop if _m == 2
	drop _m
	ren count PHOTOCOUNT
	distinct metro_id if PHOTOCOUNT!=0 & PHOTOCOUNT!=.
	local N_photo_c=r(ndistinct)
											
* MERGE POPDENS	
	merge 1:1 grid_x grid_y using "$data/125 GLOBAL CITIES/POP/GRID_POP",     
	drop if _m == 2
	drop _m		
	distinct metro_id if popdens!=0 & popdens!=.
	local N_popdens_c=r(ndistinct)
											
* MERGE COWORKING	
	merge 1:1 grid_x grid_y using "$data/125 GLOBAL CITIES/COWORKING/GRID_COWORING.dta"
	drop if _m == 2
	drop _m
	replace coworking_count = 0 if coworking_count == .
	distinct metro_id if coworking_count!=0 & coworking_count!=.
	local N_popdens_c=r(ndistinct)
											
* MERGE COUNTRY
	merge m:1 metro_id using "$data_125cities/METRO_LEVEL_COVARIATES/metro_country.dta"
	drop _m
										
* RESTRICT THE DATA TO DEVELOPABLE AREA
	drop if devle  == 0
									
* SET MISSINGS TO ZERO
	foreach var of varlist PHOTOCOUNT  {
	replace `var' = 0 if `var' == .
	}

// SET TWITTER COUNTS TO ZERO FOR METROS FOR WHICH WE HAVE TWITTER DATA	
	sum metro_id
	local max = r(max)
	forval num = 1/`max' {
		qui sum twittercount if metro_id == `num'
		local obs = r(N)
		if `obs' > 0 {
			qui replace twittercount = 0 if twittercount == . & metro_id == `num'
			display "metro `num' updated"
			}
		else {
			display "metro `num' not updated"
			}
			}
			// NOTE: TWITTER DATA ONLY FOR 70 CITIES 
									
* DROP METROS THAT WERE MANUALLY UNSELECTED
	foreach num of numlist 18 39 54 113 {
		drop if metro_id == `num'
		}	

* Save data set		
	save "$temp/gradientdata"	, replace
					
* Use data set					
	u "$temp/gradientdata", clear
								 
* Gen bin dummies and label them
	* local binsize 
	local binsize_1km=.200
	* x-limit of gradient (this is the upper bound of the interval; the bin after this interval is used as a reference bin)
	local max_x_axis=21.00
	*  starting bin 
	local drop_inner=-.4 
	g DISTB=.
	forvalues x=`drop_inner'(`binsize_1km').8{
		replace DISTB=`x' if distance_toPLborder>`x' & distance_toPLborder<=`x'+`binsize_1km'
	}
* collate all points that are below minimum into most inner bin
	replace DISTB=`drop_inner' if distance_toPLborder<=`drop_inner'
	replace DISTB=floor(distance_toPLborder) if distance_toPLborder>1
* dropping observations deep within PL for lack of sufficient observations (add 10 meters to avoid precision problems)
	drop if distance_toPLborder<`drop_inner'	
	drop if distance_toPLborder>`max_x_axis'
	tab DISTB  , gen(PLDB)
		local N_bins=r(r)
		local N_bins_less_ommited=`N_bins'-2
	tab DISTB  if distance_toPLborder<1
		local binsunder1km=r(r)
		local binsover1km=`binsunder1km'+1
* Label									
	forval bin = 1/`binsunder1km' {
		local lower =   (`bin'-1)*`binsize_1km'+`drop_inner'
		local upper = (`bin')*`binsize_1km'+`drop_inner'  
		label var PL`P'DB`bin' "`lower' < distance (km) <= `upper'"  
	}
	forval bin = `binsover1km'/`N_bins' {
		local lower =   (`bin'-`binsunder1km')
		local upper = 	(`bin'-`binsunder1km')+1  
		label var PL`P'DB`bin' "`lower' < distance (km) <= `upper'"
	}
										 
* GEN LOG COMMERCIAL DENSITY
	gen rlcomdens = log(USE_COMMERCIAL+1)-log(USE_RESIDENTIAL+1)
									
* VARIABLES OF INTEREST
	gen loc= log(any_count+1)
		label var loc "Ln SNL and Starbucks count + 1"
	gen lSNL= log(snl_count+1)
		label var lSNL "Ln SNL-S&P investment count + 1"
	gen lSTAR= log(star_count+1)
		label var lSTAR "Ln Starbucks count + 1"			
	gen ltc= log(twittercount+1)
		label var ltc "Ln Twitter count + 1"
	gen lsqft = log(snl_space)
		label var lsqft "Ln floor space"
	gen locc = log(occupancy_rate)
		label var locc "Ln occupancy rate"
	gen lrent = log(rent)
		label var lrent "Ln rent"
	gen lvalue = log(appraisal_value)
		label var lvalue "Ln appraisal value"
	gen lheight = log(height)
		label var lheight "Ln building height"
	gen lheight_activity_2000 = log(height_activity_2000)
		label var lheight_activity_2000 "log height: new constuctions"			
		label var rlcomdens "Ln no commercial - Ln no residential buildings"	
	gen lPHOTOCOUNT = log(PHOTOCOUNT+1)	
		label var lPHOTOCOUNT "Ln photo density"		
	gen DISTKM=DISTB
		label var DISTKM "Distance from prime location (km)"
	gen DISTKM_SQ = DISTKM^2
		label var DISTKM_SQ "Distance from prime location (km) squared"
	gen lcoworking_count  = ln(coworking_count+1)
		label var lcoworking_count "Ln coworking spaces"
	gen lEWPP = ln(empl+1)
		label var lEWPP "Ln emp. density"					
	label var rent "Rent"
	label var height "Building height"
	label var height_activity_2000 "Building height: Recent constructions"
	label var lpopdens "Ln population density"									
								
* Define outcomes of interest	
	global valvarlist lEWPP lcoworking_count lSNL lSTAR lpopdens ltc lPHOTOCOUNT     lheight   rlcomdens  
										
* Summary statistics for how many establishments per metro 		
	foreach x in  $valvarlist {
		egen howmany_`x'=max(`x'), by(metro_id)
		replace howmany_`x'=. if howmany_`x'==0
		di in r "`: var la `x''"
		distinct metro_id if howmany_`x'!=. 
	}
										
// DEFINE City Structure **************************************************
	gen MONO=1 if  N_PLs_city==1
	gen DUO= 1 if N_PLs_city==2
	gen POLY=1 if MONO==. & DUO==.
	foreach x in MONO DUO POLY{
		replace `x'=0 if `x'==.
	}
	global cont MONO DUO POLY
								
* RUN BIN REGRESSIONS WITH METRO FIXED EFFECT
	local N_bins_less_ommited=26
	foreach cont in $cont {
		foreach var of varlist $valvarlist {
			eststo: areg `var' PLDB1-PLDB`N_bins_less_ommited' if devle == 1 & `cont' == 1 , abs(metro_id) cluster(metro_id)
				estadd scalar metro_areas = e(N_clust)
				estadd local metro_FE = "Yes"
				forval bin = 1/`N_bins_less_ommited' {
					scalar b`var'_`bin'_`cont' = _b[PLDB`bin'] 
					scalar se`var'_`bin'_`cont' = _se[PLDB`bin'] 
				}
		}
		eststo clear							
	}
									  
* WRITE A NEW DATA SET FOR GRAPHICAL ILLUSTRATION	
	local valvarlist  lEWPP lcoworking_count lSNL lSTAR lpopdens ltc lPHOTOCOUNT  lrent  lheight lheight_activity_2000 rlcomdens   
	clear 
	set obs `N_bins_less_ommited'
	gen bin = _n
	foreach cont in $cont {
		foreach var in $valvarlist {
			gen b`var'_`cont' = .
			gen se`var'_`cont' = .
		}
	}
									
* READ SCALARS 
	foreach cont in $cont {
		foreach var in $valvarlist {
			forval bin = 1/`N_bins_less_ommited' {
				replace  b`var'_`cont' =  b`var'_`bin'_`cont' if bin == `bin' 
				replace  se`var'_`cont' =  se`var'_`bin'_`cont' if bin == `bin'
			}
		}
	}

* COMPUTE CIs
	foreach cont in $cont {
		foreach var in $valvarlist {
			gen cup`var'_`cont' = b`var'_`cont' + 1.645* se`var'_`cont'
			gen clo`var'_`cont' = b`var'_`cont' - 1.645* se`var'_`cont'
		}
	}
									
* GEN DIST
	gen dist= ((bin-1)*`binsize_1km'+`binsize_1km'/2)-(-`drop_inner') if bin<=`binsunder1km'
	replace dist=bin-`binsunder1km'+.5 if bin>`binsunder1km'
									
* Add a reference bin for visualisation
	local newN = _N+1
	set obs `newN'
	replace dist = 20 if dist == .
	foreach var of varlist cup* clo* b* {
		replace `var' = 0 if `var' == .
		}
	gen  ln_dist=log(1+dist)
	
* Generate gradient plots for outcomes of interest	
	foreach var in $valvarlist {
		if "`var'" == "loc" {
			local label "ln(SNL and Starbucks count)"
			}
			else {
			}
		if "`var'" == "lSNL" {
			local label "ln(SNL-S&P investment density)"
			}
			else {
			}
		if "`var'" == "lSTAR" {
			local label "ln(Starbucks density)"
			}
			else {
			}
		if "`var'" == "lpopdens" {
			local label "ln(population density)"
			}
			else {
			}
		if "`var'" == "ltc" {
			local label "ln(tweet density)"
			}
			else {
			}
		if "`var'" == "lsqft" {
			local label "ln(SNL office space)"
			}
			else {
			}
		if "`var'" == "locc" {
			local label "ln(SNL occupancy rate)"
			}
			else {
			}
		if "`var'" == "lrent" {
			local label "ln(office rent)"
			}
			else {
			}
		if "`var'" == "lvalue" {
			local label "ln(SNL appraisal value)"
			}
			else {
			}
		if "`var'" == "lheight" {
			local label "ln(building height)"
			}
			else {
			}
		if "`var'" == "lheight_activity_2000" {
			local label "ln(building height: new constructions)"
			}
			else {
			}
		if "`var'" == "rlcomdens" {
			local label "ln(no commercial) -  ln(no residential buildings)"
			}
			else {
			}
		if "`var'" == "lPHOTOCOUNT" {
			local label "ln(photo density)"
			}
			else {
			}
		if "`var'" == "lcoworking_count" {
			local label "ln(co-working density)"
			}
			else {
			}	
		if "`var'" == "lEWPP" {
			local label "ln(employment density)"
			}
			else {
			}			
											
		twoway /// 
				(rarea  cup`var'_MONO clo`var'_MONO ln_dist if dist<0, sort color(red%15) lwidth(none) ) ///
				(rarea  cup`var'_MONO clo`var'_MONO ln_dist if dist>0, sort color(red%15) lwidth(none)) ///
				(line b`var'_MONO ln_dist if dist<0, lpattern(shortdash) lcolor(red)) ///
				(line b`var'_MONO ln_dist  if dist>0, lpattern(shortdash) lcolor(red)) 	///
				(rarea  cup`var'_DUO clo`var'_DUO ln_dist  if dist<0, sort color(blue%15) lwidth(none)) ///
				(rarea  cup`var'_DUO clo`var'_DUO ln_dist  if dist>0, sort color(blue%15) lwidth(none) ) ///
				(line b`var'_DUO ln_dist  if dist<0, lpattern(longdash) lcolor(blue)) ///
				(line b`var'_DUO ln_dist  if dist>0, lpattern(longdash) lcolor(blue)) ///
				(rarea  cup`var'_POLY clo`var'_POLY ln_dist  if dist<0, sort color(black%15) lwidth(none) ) ///
				(rarea  cup`var'_POLY clo`var'_POLY ln_dist  if dist>0, sort color(black%15) lwidth(none) ) ///
				(line b`var'_POLY ln_dist  if dist<0, lpattern(longdash) lcolor(black) ) ///
				(line b`var'_POLY ln_dist  if dist>0, lpattern(longdash) lcolor(black) ) ///
				,  xlabel(-0.287 "-.25" 0 "0" 0.693 "1" 1.253 "2.5" 1.79 "5" 2.34 "10"  3.045 "20") ylabel(#3) graphregion(color(white) lwidth(thick))  yline(0) xline(0) /// 
				xtitle("Distance from prime location (km)") title("`label'",  size(medium)) name(`var', replace) ///
				legend(order(2 "Moncentric" 6 "Duocentric" 10 "Polycentric") cols(3) size(small)) // 
		graph save "$temp/_FIG_`var'", replace		
	}
																							
* Combine graphs
	grc1leg "$temp/_FIG_lEWPP" ///
			"$temp/_FIG_lSNL" ///
			"$temp/_FIG_lcoworking_count" ///
			"$temp/_FIG_rlcomdens" ///
			"$temp/_FIG_lSTAR" ///
			"$temp/_FIG_lheight" ///
			"$temp/_FIG_lpopdens" ///
			"$temp/_FIG_ltc" ///
			"$temp/_FIG_lPHOTOCOUNT" ///
			, graphregion(color(white) lwidth(thick)) scale(0.75)  
* Make directory and save
	capture mkdir "$figures_App/GlobalCities"
* Save Appendix Figure B.4.1.	
	graph export  "$figures_App/GlobalCities/FIG_B4_1_PLGRADIENTS_CITYSTRUCTURE.png", width(2400) height(1800) replace	
										
										
* Part II: Gradients by world region *******************************************

* Basline grid data
	u "$temp/125CitiesPLs/grid125_PL_output.dta", clear
	
* Merge distances
	merge 1:1 grid_x grid_y using "$temp/grid125_PL_distances"
	drop if _m == 2
	drop _m
	
* Merge SNL and Starbucks data
	merge 1:1 grid_x grid_y using "$data_125cities/STARBUCKS_SNL/grid_StarSNL.dta"
	drop if _m == 2
	drop  _m
	distinct metro_id if snl_count!=0 & snl_count!=.
	local N_snl_c=r(ndistinct)
	distinct metro_id if star_count!=0 & star_count!=.
	local N_star_c=r(ndistinct)
								
* MERGE TWITTER DATA
	merge 1:1 grid_x grid_y using "$data_125cities/TWITTER/DATA.dta",  
	tab _m
	drop if _m == 2
	drop _m
	distinct metro_id if twittercount!=0 & twittercount!=.
	local N_tc_c=r(ndistinct)

* MERGE EMPORIS DATA 1900 AND TODAY
	merge 1:1 grid_x grid_y using "$data_125cities/EMPORIS/EMPORIS_VALIDATION_maxHeight.dta",  
	tab _m
	drop if _m == 2
	drop _m
	distinct metro_id if twittercount!=0 & twittercount!=.
	local N_tc_c=r(ndistinct)
								
* Merge photo data
	merge 1:1 grid_x grid_y using "$data_125cities/PHOTOS/PHOTOCOUNT_250M.dta" ,   	
	tab _m
	drop if _m == 2
	drop _m
	ren count PHOTOCOUNT
	distinct metro_id if PHOTOCOUNT!=0 & PHOTOCOUNT!=.
	local N_photo_c=r(ndistinct)
								
* MERGE POPDENS	
	merge 1:1 grid_x grid_y using "$data_125cities/POP/GRID_POP",     
	drop if _m == 2
	drop _m		
	distinct metro_id if popdens!=0 & popdens!=.
	local N_popdens_c=r(ndistinct)
								
* MERGE COWORKING	
	merge 1:1 grid_x grid_y using "$data_125cities/COWORKING/GRID_COWORING.dta"
	drop if _m == 2
	drop _m
	replace coworking_count = 0 if coworking_count == .
	distinct metro_id if coworking_count!=0 & coworking_count!=.
	local N_popdens_c=r(ndistinct)
								
* MERGE COUNTRY
	merge m:1 metro_id using "$data_125cities/METRO_LEVEL_COVARIATES/metro_country.dta"
	drop _m
							
* RESTRICT THE DATA TO DEVELOPABLE AREA
	drop if devle  == 0
						
* SET MISSINGS TO ZERO
	foreach var of varlist PHOTOCOUNT  {
	replace `var' = 0 if `var' == .
	}

// SET TWITTER COUNTS TO ZERO FOR METROS FOR WHICH WE HAVE TWITTER DATA	
	sum metro_id
	local max = r(max)
	forval num = 1/`max' {
	qui sum twittercount if metro_id == `num'
		local obs = r(N)
		if `obs' > 0 {
			qui replace twittercount = 0 if twittercount == . & metro_id == `num'
			display "metro `num' updated"
		}
		else {
			display "metro `num' not updated"
		}
		}
		// NOTE: 	TWITTER DATA ONLY FOR 70 CITIES 
						
* DROP METROS THAT WERE MANUALLY UNSELECTED
	foreach num of numlist 18 39 54 113 {
	drop if metro_id == `num'
	}	
						
* Gen bin dummies and label them
	* local binsize 
	local binsize_1km=.200
	* x-limit of gradient (this is the upper bound of the interval; the bin after this interval is used as a reference bin)
	local max_x_axis=21.00
	*  starting bin 
	local drop_inner=-.4 
	g DISTB=.
	forvalues x=`drop_inner'(`binsize_1km').8{
		replace DISTB=`x' if distance_toPLborder>`x' & distance_toPLborder<=`x'+`binsize_1km'
	}
	* collate all points that are below minimum into most inner bin
	replace DISTB=`drop_inner' if distance_toPLborder<=`drop_inner'
	*drop  if distance_toPLborder<=`drop_inner' 
	replace DISTB=floor(distance_toPLborder) if distance_toPLborder>1
						
	* dropping observations deep within PL for lack of sufficient observations (add 10 meters to avoid precision problems)
	drop if distance_toPLborder<`drop_inner'
	drop if distance_toPLborder>`max_x_axis'
	
* Label	
	tab DISTB  , gen(PLDB)
	local N_bins=r(r)
	local N_bins_less_ommited=`N_bins'-2
	tab DISTB  if distance_toPLborder<1
	local binsunder1km=r(r)
	local binsover1km=`binsunder1km'+1
	forval bin = 1/`binsunder1km' {
							local lower =   (`bin'-1)*`binsize_1km'+`drop_inner'
							local upper = (`bin')*`binsize_1km'+`drop_inner'  
							label var PL`P'DB`bin' "`lower' < distance (km) <= `upper'"  
							}			
	forval bin = `binsover1km'/`N_bins' {
							local lower =   (`bin'-`binsunder1km')
							local upper = 	(`bin'-`binsunder1km')+1  
							label var PL`P'DB`bin' "`lower' < distance (km) <= `upper'"
							}

* Finalize data set for regreesions  							
							
	* GEN LOG COMMERCIAL DENSITY
		gen rlcomdens = log(USE_COMMERCIAL+1)-log(USE_RESIDENTIAL+1)				
	* VARIABLES OF INTEREST
		gen loc= log(any_count+1)
			label var loc "Ln SNL and Starbucks count + 1"
		gen lSNL= log(snl_count+1)
			label var lSNL "Ln SNL-S&P investment count + 1"
		gen lSTAR= log(star_count+1)
			label var lSTAR "Ln Starbucks count + 1"			
		gen ltc= log(twittercount+1)
			label var ltc "Ln Twitter count + 1"
		gen lsqft = log(snl_space)
			label var lsqft "Ln floor space"
		gen locc = log(occupancy_rate)
			label var locc "Ln occupancy rate"
		gen lrent = log(rent)
			label var lrent "Ln rent"
		gen lvalue = log(appraisal_value)
			label var lvalue "Ln appraisal value"
		gen lheight = log(height)
			label var lheight "Ln building height"
		gen lheight_activity_2000 = log(height_activity_2000)
			label var lheight_activity_2000 "log height: new constuctions"			
			label var rlcomdens "Ln no commercial - Ln no residential buildings"	
		gen lPHOTOCOUNT = log(PHOTOCOUNT+1)	
			label var lPHOTOCOUNT "Ln photo density"		
		gen DISTKM=DISTB
			label var DISTKM "Distance from prime location (km)"
		gen DISTKM_SQ = DISTKM^2
			label var DISTKM_SQ "Distance from prime location (km) squared"
		gen lcoworking_count  = ln(coworking_count+1)
			label var lcoworking_count "Ln coworking spaces"
		gen lEWPP = ln(empl+1)
			label var lEWPP "Ln emp. density"	
		label var rent "Rent"
		label var height "Building height"
		label var height_activity_2000 "Building height: Recent constructions"
		label var lpopdens "Ln population density"	
* Define outcomes of interest	
	global valvarlist  lEWPP lcoworking_count lSNL lSTAR lpopdens ltc lPHOTOCOUNT     lheight   rlcomdens  // lsqft locc lvalue

* Summary statistics for how many establishments per metro 		
	foreach x in  $valvarlist {
		egen howmany_`x'=max(`x'), by(metro_id)
		replace howmany_`x'=. if howmany_`x'==0
		di in r "`: var la `x''"
		distinct metro_id if howmany_`x'!=.							  
		}
							
// DEFINE NORTH AMERICAN CITIES **************************************************
	gen NA=  country == "USA" |  country == "Canada"			
	gen nNA=  NA == 0
	global cont NA nNA
						
* RUN BIN REGRESSIONS WITH METRO FIXED EFFECT
	foreach cont in $cont {
		foreach var of varlist $valvarlist {
			eststo: areg `var' PLDB1-PLDB`N_bins_less_ommited' if devle == 1 & `cont' == 1 , abs(metro_id) cluster(metro_id)
				estadd scalar metro_areas = e(N_clust)
				estadd local metro_FE = "Yes"
				forval bin = 1/`N_bins_less_ommited' {
					scalar b`var'_`bin'_`cont' = _b[PLDB`bin'] 
					scalar se`var'_`bin'_`cont' = _se[PLDB`bin'] 
				}
		}
		eststo clear							
	}
						
* WRITE NEW DATA SET FOR GRAPHICAL ILLUSTRATION	
	local valvarlist  lEWPP lcoworking_count lSNL lSTAR lpopdens ltc lPHOTOCOUNT  lrent  lheight lheight_activity_2000 rlcomdens   
	clear 
	set obs `N_bins_less_ommited'
	gen bin = _n
		foreach cont in $cont {
			foreach var in $valvarlist {
				gen b`var'_`cont' = .
				gen se`var'_`cont' = .
			}
		}
						
* READ SCALARS 
	foreach cont in $cont {
		foreach var in $valvarlist {
			forval bin = 1/`N_bins_less_ommited' {
				replace  b`var'_`cont' =  b`var'_`bin'_`cont' if bin == `bin' 
				replace  se`var'_`cont' =  se`var'_`bin'_`cont' if bin == `bin'
			}
		}
	}

* COMPUTE CIs
	foreach cont in $cont {
		foreach var in $valvarlist {
			gen cup`var'_`cont' = b`var'_`cont' + 1.96* se`var'_`cont'
			gen clo`var'_`cont' = b`var'_`cont' - 1.96* se`var'_`cont'
		}
	}
						
* Negative distance for within PL 
	gen dist= ((bin-1)*`binsize_1km'+`binsize_1km'/2)-(-`drop_inner') if bin<=`binsunder1km'
	replace dist=bin-`binsunder1km'+.5 if bin>`binsunder1km'
	gen  ln_dist=log(1+dist)

* Add a reference bin for visualisation
	local newN = _N+1
	set obs `newN'
	replace dist = 20 if dist == .
	foreach var of varlist cup* clo* b* {
		replace `var' = 0 if `var' == .
	}						

* Generate gradient plots for all variables of interest	
	foreach var in $valvarlist {
		if "`var'" == "loc" {
			local label "ln(SNL and Starbucks count)"
			}
			else {
			}
		if "`var'" == "lSNL" {
			local label "ln(SNL-S&P investment density)"
			}
			else {
			}
		if "`var'" == "lSTAR" {
			local label "ln(Starbucks density)"
			}
			else {
			}
		if "`var'" == "lpopdens" {
			local label "ln(population density)"
			}
			else {
			}
		if "`var'" == "ltc" {
			local label "ln(Twitter density)"
			}
			else {
			}
		if "`var'" == "lsqft" {
			local label "ln(SNL office space)"
			}
			else {
			}
		if "`var'" == "locc" {
			local label "ln(SNL occupancy rate)"
			}
			else {
			}
		if "`var'" == "lrent" {
			local label "ln(office rent)"
			}
			else {
			}
		if "`var'" == "lvalue" {
			local label "ln(SNL appraisal value)"
			}
			else {
			}
		if "`var'" == "lheight" {
			local label "ln(building height)"
			}
			else {
			}
		if "`var'" == "lheight_activity_2000" {
			local label "ln(building height: new constructions)"
			}
			else {
			}
		if "`var'" == "rlcomdens" {
			local label "ln(no commercial - ln no residential buildings)"
			}
			else {
			}
		if "`var'" == "lPHOTOCOUNT" {
			local label "ln(photo density)"
			}
			else {
			}
		if "`var'" == "lcoworking_count" {
			local label "ln(co-working density)"
			}
			else {
			}	
		if "`var'" == "lEWPP" {
			local label "ln(employment density)"
			}
			else {
			}			
									
								
		twoway /// 
					(rarea  cup`var'_NA clo`var'_NA ln_dist if dist<0, sort color(blue%25) ) ///
					(rarea  cup`var'_NA clo`var'_NA ln_dist if dist>0, sort color(blue%25) ) ///
					(line b`var'_NA ln_dist if dist<0, lpattern(shortdash) lcolor(blue)) 	/// 
					(line b`var'_NA ln_dist  if dist>0, lpattern(shortdash) lcolor(blue)) 	///
					(rarea  cup`var'_nNA clo`var'_nNA ln_dist  if dist<0, sort color(red%25) ) ///
					(rarea  cup`var'_nNA clo`var'_nNA ln_dist  if dist>0, sort color(red%25) ) ///
					(line b`var'_nNA ln_dist  if dist<0, lpattern(longdash) lcolor(red)) ///
					(line b`var'_nNA ln_dist  if dist>0, lpattern(longdash) lcolor(red)) ///
					,  xlabel(-0.287 "-.25" 0 "0" 0.693 "1" 1.253 "2.5" 1.79 "5" 2.34 "10"  3.045 "20") ylabel(#3) graphregion(color(white) lwidth(thick)) yline(0) xline(0) /// 
					xtitle("Distance from prime location (km)") title("`label'",  size(medium)) name(`var', replace) ///
					legend(order(4 "North America" 8 "Rest of world") cols(3) size(small)) // ,  size(small) legend(off)
		graph save "$temp/_FIG_`var'", replace		
	}
	
* Combine graphs
	grc1leg "$temp/_FIG_lEWPP" "$temp/_FIG_lSNL" ///
			"$temp/_FIG_lcoworking_count"  ///
			"$temp/_FIG_rlcomdens" ///
			"$temp/_FIG_lSTAR" ///
			"$temp/_FIG_lheight" ///
			"$temp/_FIG_lpopdens" ///
			"$temp/_FIG_ltc"  ///
			"$temp/_FIG_lPHOTOCOUNT"  ///
			, graphregion(color(white) lwidth(thick)) scale(0.75) 
								
* Make directory
	capture mkdir "$figures/GlobalCities"
* Write Figure 3
	graph export  "$figures/GlobalCities/FIG_3_1_PLGRADIENTS_NA_ROW.png", width(2400) height(1800) replace	
							
* Script ends							
							
								